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Two-dimensional XY models with resistively shunted junction (RSJ) dynamics and time depen- 
dent Ginzburg-Landau (TDGL) dynamics are simulated and it is verified that the vortex response 
is well described by the Minnhagen phenomenology for both types of dynamics. Evidence is pre- 
sented supporting that the dynamical critical exponent z in the low-temperature phase is given by 
the scaling prediction (expressed in terms of the Coulomb gas temperature t'^'^ and the vortex 
renormalization given by the dielectric constant e) z ^ - 2 > 2 both for RSJ and TDGL 

and that the nonlinear IV exponent a is given by a = z -f 1 in the low-temperature phase. The 
results are discussed and compared with the results of other recent papers and the importance of 
the boundary conditions is emphasized. 



PACS numbers: 74.50-hr, 74.40-hk, 74.25.Fy, 74.76.-w 



I. INTRODUCTION 

Superconducting films and two-dimensional (2D) 
Josephson junction arrays as well as "^He films undergo 
Kosterlitz-Thouless (KT) type transitions fromj-tiie su- 
perconducting/superfluid to the normal state.Ell3 The 
KT transition is driven by thermally created vortejfe 
antivortex pairs which start to unbind at the transition.tl 
This means that some dominant characteristic features 
of the physics close to the transition are associated with 
vortex pair fluctuations. The great current interest in 2D 
vortex fluctuations stems from the fact that they are also 
present in high-Tj, superconductors, not only in the case 
of thiit,films, but also in 3D samples just above the tran- 
sition.B It is therefore of interest to understand the prop- 
erties associated with these thermally created vortices. 
Whereas there is a fairly good consensus on the static 
properties associated with vortex pair fluctuationsja the 
dynamical aspects are less clear and some features are 
still controversial. 

The knowledge of the dynamical properties of vor- 
tex fluctuations mainly comes froni.£xperiments on su- 
perconducting fi\io& and ^He films,a0 and from various 
model simulations .B The theoretical attempts are so far 
on a pather phenomenological levelBou with few excep- 
tions The more explicit knowledge derives from several 
kinds of simulations: XY models with time dependent 
Ginzburg-Landau (TDGL) dynamics,LJ XY models with 
resistively shunted Josephson junction (RSJ) dynam=. 
icSjEliJ the Coulomb gas model with Langevin dynamics,E£l 
and the lattice Goulomb gas model with Monte Garlo 
dynamics. Lj There exist two phenomenological descrip- 
tions: the jAmbegaokar-Halperin-Nelson-Siggia (AHNSI 
descriptionci and the Minnhagen phenomenology (MP).El 
There are, likewise, two distinct proposals for the non- 
linear IV exponent a, i.e., cahns (Ref- §) and Oscaic 
(Ref. 12) with a corresponding proposal for a critical dy- 
namical exponent z = flgcaio — 1 (Ref. |l^) in the low- 
temperature phase. It has also been argued that the 
nonlinear IV exponent with the value Cscaio applies to 



an intermediate current range whereas (Mhns should be 
recovered in the true small-current limit. eI This argument 
rests on the assumption that for any finite current there 
are free vortices present and furthermore that these free 
vortices cap. be described by a conventional dynamics 
with z = 2.H 

In this paper we present extensive simulations of 2D 
XY models with RSJ as well as TDGL dynamics using 
an unconventional boundary condition. This enables us 
to obtain more information on the vortex dynamics for 
these models. 

The situation is roughly the following: The MP form 
of the dynamical response gives a goocLdescription of the 
2D XY models with TDGL dynamics ,D the Coulomb gas 



model with Langev: 
2D superconductors 



ijaiamics]i3 and experiments on 
!HU In the present paper we show 
that it also gives a good description of 2D XY mod- 
els with RSJ dynamics. The dynamical exponent z for 
the lattice Coulomb gas with Monte Carlo dynamics has 
from simulatipas been inferred to have the scaling value 
z ~ flscaic — 1. til In the present paper we verify this result 
for the XY models with both RSJ and TDGL dynam- 
ics. This is seemingly in contradiction to the results in 
Ref. I that the 2D XY models with RSJ and TDGL dy- 
namics behave differently and appear to have different z 
values. The nonlinear IV exponent a has been found to 
have the scaling value ascaic for the Coulomb gas with 
Langevin dynamical3 and the lattice Coulomb gas with 
Monte Carlo dynamics.EJ However, contradictory results 
have been found for the XY model with RSJ dynamics, 
e.g., a = OAHNS in Ref- §| and a = ascaie in Ref. In 
the present paper we find support for a = Oscaic for the 
2D XY model with RSJ dynamics. 

The picture emerging from our perspective is a generic 
vortex response well described by the MP form of the fre- 
quency response, the scaling exponent ttscaic and the cor- 
responding dynamical exponent z = Ogcaie — 1 • According 
to our view this generic vortex response describes both 
Coulomb gas models and 2D XY models and is insensi- 
tive to the detailed type of the dynamics be it Coulomb 
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gas Langevin-, Monte Carlo-, TDGL-, or RSJ-type. 

The content of the present paper is the following: In 
Sec. H we describe the XY-type models and the relevant 
correlation and response functions, as well as the relation 
to the vortex and Coulomb gas degrees of freedom. We 
also discuss the validity of linear response and the rela- 
tion between the complex impedance and the dielectric 
function of the Coulomb gas. In Sec. II] the dynami- 
cal equations are described and the boundary condition 
is introduced and discussed. Sections [V and ^ contain 
our simulation results; Sec. IV the equilibrium ones and 
Sec. ^ the result when the system is driven by an exter- 
nal current. Finally in Sec. VI we summarize our results 



and make some final remarks. 



II. XY MODEL 

On a phenomenological level, a 2D superconduc- 
tor/superfluid can be described by an order parameter 
-0(r) — |^(r)|e'^'^''', where |'(/'(r)P is proportional to the 
superfluid densitAj and V6'(r) is proportional to the su- 
perfluid velocity.B The energy associated with the order 
parameter is the kinetic energy of the current and co 
sequently the energy is proportional to / d^r[V9{r)]^ /2 
A positive (negative) vortex centered at a certain point 
is associated with the topological excitation character- 
ized by that the line integral J V6'(r) • dl of an arbitrary 
small closed loop around the point is equal to 27r (— 27r). 
There is a precise mapping between the vortices of a 2D 
superconductor and 2D Coulomb gas charges .□ Since our 
interest in the present paper is the dynamical effects of 
the thermal vortex fluctuations, we will describe our re- 
sults in the language of 2D Coulomb gas charges. 

The Xy-type models in a broad sense are models 
representing the continuum order parameter ipir) = 
|V'(r)|e*'''-'''' put on a lattice. Let us for convenience 
choose a square lattice. The discretized version is then 
ijjj = IV'jle**^ , where the index j denotes the lattice 
points. Let us simplify further by neglecting the vari- 
ations of the magnitude of the order parameter and take 
= to be a constant. The discretized version of 
the energy then takes the form 



H 



XY 



= J 



U{4>^j = 0, -9,), 



(1) 



where J < 



is termed the XY coupling constant and 



the sum is over nearest-neighbor pairs. The lattice con- 
stant is taken to be unity so that (pij = 6i — 9j corresponds 
to W6 (in the direction from j to i). The function C/(0) 
has to be equal to 0^/2 for small (j) in order to yield the 
correct continuum limit and in addition U{(f)) has to be a 
periodic function of 27r since the phase angle 9i for each 
lattice point is only defined upto a multiple of 27r. A 
possible choice for U{(j)) is then 

U{(j)) = 1 - cos^ 



and with this choice the model is the usual 2D XY model 
or the planar rotor model. This particular interaction 
would, e.g., arise if each lattice point was a small su- 
perconducting island which was Josephson coupled to its 
nearest neighbors, and the system is called a Josephson 
junction array (JJA). We will use this choice of the in- 
teraction in the present paper. However, from the point 
of view of vortex fluctuations any U{4>) fulfilling the nec- 
essary requirements stipulated above is a valid choice. A 
possible generalization is 



cos 



2p' 



(2) 



where p = 1 corresponds to the usual XY model. The 
practical point with such a generalizatiop-js that the vor- 
tex density increases with increasing pS3 Consequently 
the vortex response is sometimes easies to extract from 
simulations for a p value larger than 1 J3 

The Boltzmann factor for a particular configuration is 
given by e"^-^'^/-^ where T is the temperature in units 
of fcs = 1. From this all thermodynamic properties can 
be obtained. 

The mapping between the XK^odel and the Coulomb 
gas representation is as follows:li3 The effective temper- 
ature variable for the Coulomb gas charges is given by 
yCG ^ T/[2nJ{U")], where T is the temperature for 
the XY model, (• • •) denotes a thermal average, and 
U" = d'^U/dcjP. The supercurrent through a link is given 
by JU' = JdU /d(j>. The Coulomb gas charge ri/ , corre- 
sponding to an elementary plaquette of the square lattice 
I, is given by the directed sum (corresponding to a line 
intsgral) over the four links (ij) making up the plaque- 
tte:N 



ni 



T 



The correlation fimction G{k,t) is a key quantity and is 
defined by 

where F{k,t) is the ID Fourier transform 
F(fc,t)=^F,„(t)e^'=™, 

m 

m labels the rows of the lattice, and finally 

{ij}em 

where the summation is over all the links making up the 
row m. The Fourier transformation of the charge density 
correlation function g{k,t) is related to G{k,t) by 
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f T Yg{k,t) 



(3) 



Linear-response theory then hnks olk^t) with the dielec- 
tric response function l/e{k,uf) byQ 



Re 



t{k,u) 



2nLuT 



CG (-oo 



i{k,o) 



Im 

where 



i{k,Lu)_ 



2t:luT 



rp2 



CG /-oo 



dt sinLut G{k, t), 
(4) 



y2 



dt COS iLjtG{k,t), (5) 



' l-^G(fc,0). 



i{k,o) 



2^2 



(6) 



The quantities l/e{0,uj) and (5(0, i) will be of particular 
interest in the present investigation. 

The thermodynamic KT transition is characterized by 

1 1 

fe^o e(k, 0) e 



below the transition and 



lim 



1 



k^o i{k,0) 







above. Precisely at the transition lim^^o V^(^7 0)r'nr-i 
jumps from the universal value l/eT^*^ = 4 to zeroO'tS 
The equal-time correlations fall off like power laws with, 
distance below the transition and exponentially above.0 
For example, the correlation function G{r,t — 0) falls off 
like 



(7) 



below the transition temperature. The fact that the cor- 
relations decay algebraically with distance reflects that 
the whole low-temperature phase is quasicritical. 

As explained in the previous section one motivation for 
the present paper is the question of the generality of t 
MP form for the dynamical response, which is given b 



KT transition whereas it can only be approximately cor- 
rect above because of the presence of free vortices which 
always dominates the response for small enough fcequen- 
cies and gives a Drude-like response in this limit .0 In the 
present paper we focus on the low-temperature phase. 
In this case the leading small uj behavior of Eqs. (H) 
and (H) refl£cts a 1/t decay for large t of the function 
may also observe that Eq. (^) leads to 
a logarithmic divergence of the real part of the conduc- 
tivity: (7(w) — Ci;Im[l/e(fc = 0,a;)] —Inuj for small 
which is compatible with standard scaling argument by 
Fisher and Fisher, Fisher, and Huse in Ref. [l9ll^ 

The two features G{r,t = 0) oc r-I^i/""^ ^nd 
G{k = O^t) — J d^rG{r,t) oc 1/t can be turned into 
an argument for the dynamical critical index z in the 
following way£3 We assume that G'(r, t) must be of the 
form 

G{r,t)<xX"fir/X,t/T,a/r,Ta/t), 

where A is the correlation length or screening length 
which diverges in the low-temperature phase, r is the 
corresponding diverging relaxation time so that 

T oc A^, 

where z is the dynamical exponent. In addition we have 
a short distance scale a, i.e., the lattice constant or the 
size of a Coulomb gas particle and a nondiverging charac- 
teristic time scale r^, i.e., Tq oc P/D where Z? is a vortex 
or Coulomb particle diffusion constant and I is some non- 
diverging length scale like ^ = a or ^ = 1 / where n is 
the density of Coulomb gas particles. Let us choose t = 
and r = A so that 

G(r,0) cx r"/(l,0,a/r,cx)) 

and make the ad hoc scaling assumption that 

lim /(I, 0, a/r, oo) — /(I, 0, 0, oo) = const, 

r — >oc 

where const 7^ and 7^ ±00. This requires a = 
_^ 2 since G(r,0) cx r--[(i/?T^'=^)-2] . We then 
also have that 



d2rG(r,i) = A-(iA^^"")+^ / d\fir/X,t/r,a/r,Tjt) 



Re 



Im 



1 



1 



e{k = 0,uj) e(0,0) 



1 



UJ 



1 



e(fc = 0,uj) 



2 ujujq In uj/ujq 



en 



(8) 



(9) 



The characteristic frequency ujq vanishes as the KT tran- 
sition is approached from above and below. Q The idea 
behind the MP form is that it describes the response due 
to the bound pairs. Consequently, it is expected to have 
the correct leading small-frequency behavior below the 



Now we choose \ = t^ so that 



and assume that 

lim /(r/^l/^ 1, a/r, rji) = /(0, 1, a/r, 0) = /(a/r), 

I— »oo 

where j{x) is a well-behaved function so that 

'dVG(r,t)cxt[-(i/^^^'''''+2]/. rf2^^~(^/^) 
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for large t. This is consistent with J (PrG{r,t) cx 1/t 
provided 



- 2. 



(10) 



The dynamical exponent z given by Eq. ( p^ ) has been 
inferred through simulations pof the lattice Coulomb gas 
with Monte Carlo dynamics O In the present paper we 
conclude that the same is true for the XY models both 
with RSJ and TDGL dynamics. |_. 

It has been argued by DorseyjHil using scaling anal- 
ysis, that for a 2D superconductor the exponent a in 
the nonlinear IV characteristics V (x I'^ has the value 
a = z + 1 precisely at the KT tMnsition. It has fur- 
ther been suggested by MinnhagentS that since the whole 
low-temperature phase is quasicritical the same rela- 
tion should apply throughout the low-temperature phase. 
This together with Eq. leads to the prediction 



Z+1 



1 



- 1. 



(11) 



The nonlinear IV exponent a — flscaic in Eq. ( |ll| ) 
has been inferred through simulations-.for the Coulomb 
gas model with Langevin dynamicsE2l and theJ^ittice 
Coulomb gas model with Monte Carlo dynamics. til 

The response to an imposed current is for a 2D|SjLicer- 
conductor given by the complex impedance Z(tt'):aE3 

where E(w) is the frequency dependent electric field 
and j(w) is the current density. Or equivalently for a 
quadratic sample V{ijj) = Z{uj)I{llj), where V is the volt- 
age across the superconductor in some direction and / 
is the total current in the same direction. The linear- 
response function Z~^{oj) is related to the Coulomb gas 
linear-response function l/e{k = 0,uj) by 



Po 



iu}e{k = 0, w) ' 



(12) 



where po is the density of superconducting electrons 
which for an XY model is given by J{U"). This means 
that the effect on the vortex fluctuations of an imposed 
current is given by l/e(fc = 0,lo). For small oj this is the 
dominant contribution. 

It is instructive to consider the linear response to an 
imposed current directly in the the case of the XY model 
with RSJ dynamics. Let us consider a quadratic lattice 
and let {ij)x be a link at position r parallel to the x axis 
and denote the difference in phase angle by (f)ij — WxO{r); 
when the coupling to the electromagnetic field is included 
(jjij denotes the gauge invariant phase difference. The su- 
percurrent through the link at time t is JU'[\7 ^9{r,t)] 
and the normal current is proportional to —\7a:0{r,t) 
where the dot denotes the time derivative. Thus the total 
current ix{j^, t) through the link is 



ix{r,t) = -V,^(r,i) + JU'[Vxe{r,t)] 



(13) 



in some convenient unit system. The voltage in the RSJ 
model is proportional to the normal current so we can de- 
fine the response function corresponding to the complex 
impedance as Z{r — r' ,t — t') = P{r — r' ,t — t'), where 



P{r-r',t-t') = 



d{Vx0{r,t)) 



dix{r',t') 



(14) 



It is shown in the appendix that the Fourier transform 
of P is given by 



PCk,Lu) 



where po = J{U") so that 



Po 



e(k,^) 



Z(k,w) = 



1 Po 
iuj e(k, Lo) 



(15) 



(16) 



This means that the response to a uniform time vary- 
ing current is given by Z{uo) — ^(0,w). Below the KT 
transition we have 

lim lim — = oo 

w^ok^o e(k, Lj) 

so that the static response to a uniform static current 
below the KT transition is nonlinear. However, for 
any finite frequency the response is linear to the low- 
est order. One also notes that in the limit of high fre- 
quency l/ia;e(k, w) vanishes and Z in Eq. ( p^ reduces 
to Z{'X)) — 1, which means that the response in this 
limit is given by the resistive shunt in the RSJ model. 
For smaller frequencies the response is given by the vor- 
tex fluctuation Z{lj) oc iuje[Q,uj) / as already stated in 
Eq. (|l|). 



III. DYNAMICAL EQUATIONS AND 
BOUNDARY CONDITIONS 

Simulations by necessity involve lattices with a finite 
linear dimension L from which the results for the ther- 
modynamic limit L oo have to be extracted. This 
means thatria practice the choice of boundary condition 
is essential.c3 The most commonly used boundary con- 
dition in order to extract the thermodynamic limit for 
the XY models is periodic boundary conditions (PBC) 
imposed on the phase angles 6i. However, as discussed 
in Ref. ^ the PBC for the phase angles leads to a 
nonperiodic boundary condition for the vortex interac- 
tion. The boundary condition for the phase angles which 
corresponds to a periodic vortex interaction is ipstead 
the fluctuating twist boundary condition (FTBC).E3 The 
dynamics we are investigating in the present paper are 
linked to the vortex fluctuations and consequently the 
natural boundary condition is PBC for the vortices. This 
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is the commonly used boundary condition for simulatioaa 
of the lattice Coulomb gas with Monte Carlo dynamicsEll 
andJjhe continuum Coulomb gas with Langevin dynam- 
ics.E3 Thus the important point in the present context 
is that PBC for the vortices means FTBC for the phase 
angles. The FTBC for the phase angles has so far^been 
used in connection with Monte Carlo simulations.ll3 In 
the present paper we extend the use of these boundary 
conditions to XY models with RSJ and TDGL dynam- 
ics.cj Of course the boundary condition should not mat- 
ter in the limit L ~* oo. However, we in the present 
paper find that by using FTBC for the phase angles we 
are able to extract more information from our finite L 
simulations. 

In this section, we briefly review the dynamical equa- 
tions of motion for RSJ and TDGL in the case of PBC 
for the phase angles. Then we construct the equations 
of motion for FTBC starting from total current conser- 
vation and the condition that the equations of motion 
should lead to the correct equilibrium distribution. We 
focus on the ordinary XY model, which corresponds to 
p = I case in the previous section, but the extension to 
a general p is straightforward. 

We begin with an LxL array of the resistively shunted 
junctions with PBC in both directions. In the RSJ dy- 
namics of 2D XY model the net current from site i to site 
j is written as the sum of the supercurrent, the normal 
resistive current, and the thermal noise current: 



ir sm 



— + F- 

r 



where ic = 2eJ/h is the critical current of the single 
junction, Vij is the potential difference across the junc- 
tion, r is the shunt resistance, and the phase angles 6i 
are periodic in both directions (6'j = Oi^Lx — Oi^^y). 
The thermal noise current Fy at temperature T is re- 
quired to satisfy (Fy(i)) = and (Fy(t)Ffci(O)) = 
(2kBT/r)6{t){SikSji — SuSjk)- The current-conservation 
law at each site, together with the Josephson relation 
d{9i — Oj)/dt = 2eVij/h, allows us to write the equations 
of motion in the form 



smi 



(17) 



where the primed summation is over four nearest neigh- 
bors of j, Gij is the lattice Green function on the square 
lattice with PBC, r]jk is the dimensionless thermal noise 
current defined by rijk = Tjk/ic, and the unit of time is 
h/2eric. The thermal noise current satisfies {i]ij{t)) — 
and 



iVij it)mi(>)) = 2T{5ik5ji - SiiSjk)S{t), 



(18) 



where T is in units of J/ks- 

In the TDGL dynamics with PBCj-en the other hand, 
the equations of motion are given byca 



where F is a dimensionless constant which determines the 
time scale of relaxation, H = —JJ2{ij) cos{9i — dj) is the 
Hamiltonian of the usual XY model, and 9i is periodic in 
both directions. The thermal noise term Fi(t) is assumed 
to satisfy {Ti{t)) = and (Fi(t)Fj(O)) = 2hrkBTSijS{t). 
After rescaling the time and the temperature in units of 
Ti/TJ and J /kB, respectively, the equations of motion for 
TDGL dynamics are written as 



sm 



(19) 



Ti/TJ satisfies 



(20) 



where the thermal noise term rji = 
(77, (i)) = and 

(ry,(i)r;,(0)) =2r%<5(i). 

In numerical simulations for PBC, we use Eqs. (|l^ and 
(|l9| ) for RSJ and TDGL dynamics, respectively, with 
the corresponding thermal noises satisfying Eqs. (|l8|) and 

®- . 

Next we consider the fluctuating twist boundary con- 
dition FTBC. In this case a variable A = {Ax^^y) is 
introduced andj-the phase difference (pij on the bond [i, j) 
is changed intal3 



- - i-y • A, 



(21) 



where = Vj — is a unit vector from site i to j, and 
the phase angles are periodic: 9i = 9i+Lx = 9i^Ly. In 
the study of equilibrium behaviors for FTBC using MC 
simulations— it is sufficient to know the Hamiltonian of 
the systemtj 



H 



(22) 



In dynamical simulations, on the other hand, we must 
also have equations of motion for the new variables A^, 
and Ay in addition to the equations of motion for phase 
variables 9i. 

The physical situation we have in mind is a sample 
where no current passes through the boundary. For 
the RSJ model, which has local current conservation, 
this implies the total current conservation condition 
/(ir^i(r,t) = 0, where \{v.t) — [ix{v,t)^iy{Y,t)] is the 
total current density at point r and the integral is over 
the whole sample. This condition can also be expressed 
as 



r 



= "I E 



- 9j - A^) - 77A, 



(23) 



(and the similar equation for the y direction), where the 
summation over all nearest- neighboring pairs 

in X direction, Vx is the voltage drop over the sample, 



5 



and r]^^ denotes the thermal noise current. This follows 
because the left-hand side is recognized as the normal 
current whereas the right-hand side is the negative of 
the sum of the supercurrent and the noise current. As 
discussed in connection with Eq. ( |l3| ) the voltage is by 
the Josephson relation proportional to V0(r,t). For the 
voltage across the sample this means that [see Eq. 



A, 



2e 



14 



(24) 



because the phase angles are by construction subject to 
periodic boundary conditions. Thus from Eqs. (^3|) and 
( p^ ) , we obtain the equations of motion for the twist vari- 
ables: 



dAy 
~~dt 



12 



12 ■ 



- - A, 



(25) 
(26) 



where we have again written t in units of h/2eric. Next 
a noise correlation consistent with the equilibrium con- 
dition has to be found. To this end we make the ansatz 
of a standard white-noise correlation (?yA^ (i)^Ai, (0)) — 
iv Ay (t)r] Ay (0)) — a'j^S(t) and determine the appropriate 
a'j^ in the following way: The equations of motion for the 
phase variables in FTBC are written as 



with 



and 



h, = -Y^ Gy ^ sin(0j - Ok - A.jk) 



{Vi]i't)Vki{0)) = cr'^iSikSji - SiiSjk)S{t), 



(27) 



(28) 



where cr^ = 2T [sec Eq. (p^J. From the full equations 
of motion for RSJ model in FTBC [Em- (H) - we 
arrive at the Fokker-Planck equation:E3 



dW 



d 



d 



d 



(hyW) 



i,3 



1 2 f d'^W d'^W 



2*^^ V dAl 



dAl 



where W — W{{9i},Ax,Ay]t) is the probability distri- 
bution function and 



-y 

12 



sin(6'j - 6*^ - A^), 



and the similar equation for hy. The stationary solu- 
tion, which satisfies dW/dt = 0, is of the correct form 



W = e with the Hamiltonian given by Eq. ( p^ ) pro- 
vided 



(3(7 



1 

l2' 



(29) 
(30) 



and consequently a\ = 2T/L^. 

The equation of motion for the twist variables are 
hence of the Langevin form 



(31) 



with Fa = and (?7a. (i)r,A. (0)) - (?7a„ (t)?7A„ (0)) = 
{2T/L^)5{t). 

In the TDGL model the total current conservation con- 
dition can still be imposed whereas the local current con- 
servation condition is relaxed. Thus Eqs. ( p5|) and (|2^ 
remain unaltered whereas the equations for the phase an- 
gles are simplified to [compare Eqs. (^ and (pi|)] 



(32) 



where we have used the dimensionless time t by intro- 
ducing the time unit of h/TJ as in Eq. ( M) . Just as 
for the RSJ case one finds that Fa = l/Lr and that 
the noise correlation (??a^ (t)?7A, (0)) = (?7a„ (i)?7A„ (0)) = 
(2T/i^)(5(i) leads to the correct equilibrium. To some 
extent the TDGL dynamics may be viewed as a simpli- 
fied version of the RSJ dynamics where the total current 
conservation is kept but the local current conservation is 
relaxed. Thus from this point of view it is perhaps not 
surprising that the two models (as we will see) have the 
same generic vortex dynamics. 

The twist variable A plays an important role in our 
analysis of the vortex dynamics and there exists a rather 
direct connection between the twist and the vortices: The 
electric field E(i) due to the vortex current density is 
perpendicular andJs, as a consequence of the Josephson 
relation, given bytll 

The connection between (jv{t)) and A is discussed in 
Ref. ^ when a vortex moves across the sample then the 
twist variable changes by 2tt/L. In other words, if the 
time to is associated with the movement of a vortex across 
the sample, then we get A = 27r/_Lio = 27r(v)/L^ where 
I (v) I = L/to is the vortex velocity. If there are iV„ moving 
vortices, then we obtain A — 2ti{Ni, / L'^My) = 27r(jt,) 
which leads to the relation given by Eq. 
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where Vx is the voltage drop across the sample in x di- 
rection (we obtain the similar equation for Ay). 

So far we have considered the situation when the total 
current in the sample is zero, which corresponds to no 
current passing over the boundary. Let us now consider 
the case when the total current is a constant dc current Id 
in the x direction. By following the steps from Eq. ( |2^ ) 
to Eqs. (25) and ( |2^ ) one obtains the modified equations 
of motion for the twist variable A 



dAx 1 



dt L2 



sin{9, -9j - Ax) + r/A, - id (33) 



^ = ^11 «"^(^* - - ^y) + ^A. (34) 

with id — Id/ L in units of ic- The voltage drop in the x 
direction [sec Eq. (|2j)] is given by 



Vx 



-LA, 



(35) 



with Vx in units of ric for RSJ and in units of F J/2e for 
TDGL, respectively. Thus the equations of motion in the 
presence of an externally imposed dc current Id in the x 
direction are given by Eqs. (27), (|33|), and ( |3^ ) for RSJ 
and by Eqs. (||), (||), and (||) for TDGL. 

An alternative and commonly used method in numer- 
ical simulations of the current-driven XY model is to 
impose uniform currents through the boundary in one 
direction. This requires an open boundary condition for 
the phase angles in the direction of the applied current 
and the periodic boundary condition can only be kept 
in the perpendicular directionHj This means that an 
open boundary is explicitly introduced. One advantage 
with the present method is that the periodic boundary 
conditions on the phase angles are kept and no explicit 
boundary is introduced. In the following two sections we 
present the results obtained from the dynamical equa- 
tions described in the present section both for the PBC 
and the FTBC. 



IV. SIMULATION RESULTS 

In this section, we present simulation results for the 
TDGL and RSJ dynamics with periodic boundary con- 
ditions PBC and the fluctuating twist boundary condi- 
tions FTBC. For PBC, we use Eqs. (P) and (|Ta) in the 
RSJ case and in the TDGL case Eqs 
FTBC, we use Eqs. (H) - 
(|^), and (H) for TDGL. 



For 



(|19|) and (gO[ 
for RSJ, and Eqs. (H) 



We integrate the equations of motion by discretizing 
time into small steps At. At each step the appropri- 
ate random noise, generated from a uniform distribu- 
tion, is introduced with {riij{t)'^) — 2T/ At for RSJ and 
imit)^) = 2T/At for TDGL [see Eqs. and^]. We 
want to integrate to as long times as possible.E30n the 
other hand the larger At we choose the larger is the error 



introduced by the discretization. In order to get a han- 
dle of the choice for At we use the following identity: Let 
us introduce a local variable on one particular site k. 
The Hamiltonian of the system is then 



H 



with Qi 7^ for i — k and ai — otherwise, and the 
partition function is given by 



with the inverse temperature f3 = 1/T. After a simple 
change of variable 9i + ai —^ 9i, we find that Z does in 
fact not depend on and thus 



\nZ 



- 0, 



from which we conclude that 



E u'{9k~e,) 



and thus 



T 



(36) 



provided we have defined T by the local correlations 



T 



U'{9k - > 



4([/") 



(37) 



where the summation is over four nearest neighbors (de- 
noted by j) of site k. The point is now that for a finite 
At one finds that T > T . In the present simulations we 
use the time step At = 0.01 for TDGL and At = 0.05 
for RSJ. These choices make T differ from T by less than 
3%. 

The fact that T > T for a finite time step suggests 
that the effect of the finite time step to some extent is 
equivalent to an increased temperature. We have tried to 
take this into account when analyzing quantities related, 
to 1/e by noting that for FTBC one has l/e(0,0) = 0,EI 
which means that [compare Eq. (||)] 

G((M3) 
([/") • 

Thus we can estimate an effective temperature by T'^^ = 
G(0,0)/(t/"). For example, for T = 0.80 and the time 
step At = 0.05 we for RSJ get T''^ « 0.82 whereas we 
for TDGL and the time step At = 0.01 get T""^ « 0.803. 
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A. Dynamical response functions with periodic 
boundary conditions 



G 



MP 



T2 



1 



We will first consider the vortex dynamics as reflected 
in the complex dielectric function given by Eqs. (Q) and 
(^. It has so far been established that the MP form 
Eqs. (||) and (^L-gives a good representation of the ex- 
perimental data,E3 as well as the simulation data for the 
TDGL dynamics of the XY model on a square lattice 
with p — 2 and on the triangwlar lattice with p — l,u and 
the 2D Coulomb gas model.E2l In the present investigation 
we find that the same is true for the XY model with RS J 
dynamics. This is illustrated in Fig. |^ which shows the 
real and imaginary parts of l/e(k = 0, w) — l/e(0, 0) with 
RSJ dynamics for T = 0.85. The full line in the figure 
has been obtained from a least-square fit to the MP form 
of the real part in Eq. (||) with two free parameters (e and 
ujq), and the broken line has been obtained by using the 
same values of the parameters in Eq. (the frequency 
range in Fig. ^ corresponds to 0.08 < uj/ujo < 4.7). The 
MP form has the characteristic feature that the ratio is 
|Im[l/e(0,w)]|/Re[l/e(0,a;) - l/e(0,0)] = 2/7r at the fre- 
quency where the imaginary part has its maximum. One 
sees directly in Fig. |l] (i.e., without any curve fitting) that 
the dotted vertical line is close to this maximum and it 
is hence easy to verify that the ratio is indeed close to 
2/tt. In short, our present simulations of the complex 
dielectric function confirm that the RSJ dynamics is well 
described by the MP form at temperatures below as well 
as somewhat above the critical temperature in agreement 
with what was found earlier for the TDGL dynamics in 
Ref. @.EI 

As pointed out in connection with Eqs. (^ and (H), 
the leading small u! dependence of the MP form 



Re 



1 



e(0,c^) e(0,0) 



OC U! 



and 



Im 



(X Lolnuj 



reflects that (^(k = 0, i) oc 1/t for large t. More precisely, 
since 



G(0,t) 



smojt 



Re 



/o <^ 
we find for the MP form 

J.2 



1 



1 



e(0,u;) e(0,0) 



duj, 



G^^^^(0,t) 



[Ci{uJot) sin Wot — si{uJot) cosujQt] , 

(38) 



where the cosine and the sine integrals are defined by 
Ci(a:) = — dt cost/t and si(a;) = — dtsint/t, re- 
spectively. In the limit of wot —> oo, Eq. ( p8[ ) reduces 
to 



This 1/t tail in the vortex correlations has been verified 
in Ref. |l^ for TDGL dynamics and in Ref. |l^ for the 
Coulomb gas model. We will here verify the same result 
for the RSJ dynamics. 

By necessity, the finite lattice sizes used in the sim- 
ulations introduce a finite relaxation time tq at large t 
for the zero-A; mode. By studying the lattice size de- 
pendence of G{0,t) we have found that this finite size 
induced relaxation changes the large-t decay from 1/t to 
(1/t) exp(— t/rc). In fact we have found that G(0,t) for 
finite lattices to a good approximation is of a modified- 
MP form (MMP): 



G 



MMP _ /.MP 



G^^^^exp(-t/TG). 



(39) 



Figure g shows ln[tG(0,t)] as a function of time for the 
system sizes L = 6, 8, 10, 12, 16, and 64 in case of (a) RSJ 
and (b) TDGL dynamics at T = 0.85. The full drawn 
curves are least-square fits to Eq. (p9|). As is apparent 
from Fig. ||, tG approaches a constant for large lattice 
sizes verifying that G indeed goes as 1 /t for large t both 
for RSJ and TDGL dynamics. 

The fits to the MMP form (full drawn curves in Fig. 0) 
show that lntG(0, t) goes as —t/rQ for large t. In Fig. H 
we have plotted tq [determined by the fit to Eq. (|3S|)] 
as a function of lattice size L in a log-log scale. From 
finite-size scaling we expect that in the low-temperature 
phase TG diverges as tq oc for large L where z is the 
dynamical critical exponent. This behavior corresponds 
to straight lines in Fig. ^ and the full straight lines in 
the figure suggest that the asymptotic scaling is reached 
already for relatively small L. Assuming that this is the 
case, we find from the slopes of the lines that for T = 0.85 
z w 1.6 in case of RSJ and z w 2 for TDGL. Thus the 
z values in case of PBC are different for the RSJ and 
the TDGL dynamics. This difference between RSJ and 
TDGL in case of periodic boundary conditions was also 
found by Tiesinga et al. in Ref. H, where in the temper- 
ature interval T e [1.1, 1.3] z w 2 for TDGL and z « 0.9 
for RSJ; the authors concluded that the TDGL some- 
what unexpectedly describes the osperiments on Joseph- 
son junction arrays by Shaw et alxB better than the RSJ 
model. The conclusion we arrive at is different since we 
find that for FTBC the equivalence between RSJ and 
TDGL is restored. The apparent difference in case of 
PBC appears to be a boundary effect E2l We believe that 
the physical situation in Ref. ^ and most other com- 
mon experimental situations are in fact better described 
by the FTBC. Of course, for large enough system sizes, 
intensive physical quantities do not depend on the ex- 
plicit choice of boundary condition. But the point here 
is that, because the relaxation of the zero-fe mode is de- 
scribed by a relaxation time tg which diverges for infinite 
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systems, the exponent which describes how this diver- 
gence is approached, awaears to be sensitive to the choice 
of boundary conditionBil 

We also note that for T = 0.90 we find z k, 1.6 in 
case of RSJ with PBC. This suggests that z for PBC ap- 
proaches a value less than 2 as the KT transition is ap- 
proached from below, although the numerical accuracy 
may be insufficient to make a firm conclusion. 



B. Dynamics for the fluctuating twist boundary 
conditions 

In case of FTBC the static dielectric function func- 
tion l/e(k, 0) is identically zero for k = 0, whereas 
limk^o l/e(k, 0) ^ below the KT transitionES In 
Ref. |l^ it was shown that for the Coulomb gas model 
with Langevin dynamics the function l/e(k,w) for small 
k is to good approximation giv en by the MP form. Since, 
as explained above in Sec. Ill, PBC for the vortices (as 
in Ref. [lO| ) corresponds to FTBC for the XY model we 
also expect to find the MP form for small k in the present 
case. This is illustrated in Fig. ^ which shows the real 
and imaginary parts of l/e(k, for k = (0, 2-k/L) with 
L = 64 for the XY model with RSJ dynamics. The full 
drawn and broken curves represent the MP form just as 
in Fig. |l| and the dotted line shows that the peak ratio 
is close to 2/7r. Figure || demonstrates that the imagi- 
nary part depends very little on the k value whereas the 
real part increases with decreasing k for fixed frequency. 
This behavior has also been found for the Coulomb gas 
model with Langevin dynamics (compare Figs. 11 and 12 
m Ref.|l^). Figure ^ shows how the relaxation time tq of 
(^(k, t) depends on A:: G cx (e~*/'^'^)/i for large t and tq 
diverges as k is decreased. In Ref. |l^ it was found that 
Tq (X k~^ for the Coulomb gas model with Langevin dy- 
namics. Our present convergence is not good enough for 
establishing this result, but Fig. ||(b) suggests that such 
a behavior is also consistent with the present simulations 
of the XY model with RSJ dynamics. 

Next we turn to the diverging relaxation time r and the 
dynamical critical exponent z for the case of FTBC. We 
will use the fact that in the low-temperature phase the. 
resistance i? of a finite system is proportional to 1/t.l3 
This follows because of the Nyquist formulaS 



R 



1 



2kBT 



dt{v{t)vm 



(40) 



which relates the resistance to the voltage fluctuations 
over the sample and the fact that V oc {d/dt)A(t) where 
A(p is the phase difference over the sample. Since Acj) is 
dimensionless it follows that i?jscales like l/r where r is 
the relevant characteristic time£3 In the low-temperature 
phase R vanishes in the limit of large system sizes since r 
diverges. Consequently the finite-size scaling i? cx l/r oc 
defines the value of the dynamical critical exponent 
z in the low-temperature phase. For the XY model with 



FTBC the phase difference over the sample in one direc- 
tion (let us choose the x direction) is given by A0 = LA^- 
It follows that R can be expressed as 



R 



2TQ 



[A,(e)-A,(0)]2 



(41) 



where T is in units of J /ks^ Ax{t) is the twist variable in 
the X direction at time t, and R is in units of the shunt 
resistance r of a single Josephson junction for the RSJ 
model and VJ /2eic for TDGL model, respectively. Since 
Eqs. (40) and (O) are identical in the limit of large 0, 
i.e., for 8 ^ r]^ we for practical reasons use Eq. (|4l| ) 
in the present simulations (we have used = 16000 and 
O ^ r). Figure ^a) shows the results for the XY model 
with RSJ dynamics for T =0.90, 0.85, and 0.80. The 
data are plotted as logi? against logL and to good ap- 
proximation fall on a straight line, whose slope gives an 
estimate of the critical exponent z, and we obtain z =2.0, 
2.7, and 3.3, respectively. Figure |^(b) shows the same 
features for the XY model with TDGL dynamics at the 
same three temperatures T =0.90, 0.85, and 0.80 and the 
estimated values of z « 2.1, 2.8, and 3.3 are close to the 
ones obtained for the RSJ dynamics. 

Thus for the FTBC we find the same z below the KT 
transition for RSJ and TDGL dynamics, which is in con- 
trast to the PBC case where we found different values of 
for each dynamics (compare the discussion of Fig. ^ in 



Sec. rV). Furthermore for FTBC we find that z appar- 
ently approaches 2 when the KT transition is approached 
from below (T = 0.90 is very close to the KT transition 
temperature) for both dynamics; this did not seem to 
be true for the RSJ dynamics with PBC [z « 1.6 at 
T = 0.90). Our conclusion is that the dynamical critical 
exponent z is a boundary sensitive quantity. We also note 
that the FTBC is adequate for describing an open sys- 
tem with voltage fluctuation across the system and that 
consequently the z values obtained for this case describe 
the most usual physical situation. 

It is in fact possible to estimate the characteristic time 
T very directly since the variable A^^ changes by the 
amount 27r/i when a vortex moves across the system 
in the y direction, as discussed in Sec. [11. Every such 



event consequently is signaled by a 27r step in the time se- 
ries of the variable LA^. Figure [s] illustrates this for the 
RSJ dynamics at T = 0.85 for various system sizes. As 
seen in the flgure the 2it jumps are very well observable. 
The characteristic time scale r of these 2it jumps is eas- 
ily estimated as the average time between the jumps and 
we expect that t ^ with the same dynamical critical 
exponent z as in i? ^ L^^ . Figure ^(a) shows r plotted 
against the system size L in a log-log plot for the RSJ 
dynamics for three different temperatures (in practice we 
use a coarse graining of 100 time units in our estimate of 
the average time between the 27r jumps). The full drawn 
straight lines in Fig. ^(a) have the slopes given by the z 
values determined previously from the calculation of R 
[see Fig. 0(a)]. As seen the two ways of determining z 
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agree very well. Figures 0(b) and |(b) illustrate the same 
agreement in case of TDGL dynamics. 

Let us now consider what happens when a finite cur- 
rent is applied across the system. The scaling argument 
by Dorseyc^ makes use of the fact that the pirrent den- 
sity J introduces a new length scale l/J.Ea This new 
length scale replaces thej-£nite size L in the leading L 
dependence of i?, so thatcJ 



1 



/oc J^J 



and consequently V oc below the KT transition 

as suggested in Rcf. [ij. From the finite-size scaling 
of R and r we obtain z and using the scaling argu- 
ment this z is related to the nonlinear lY exponent by 
a — z + 1 where 1/ - J'^. In Table | we have given the 
values of a = 2 -f 1, where z has been obtained from 
the finite-size scaling of R. Another scaling argumenttj 
gives [see Eq. (M)] z = — 2 and consequently 

a = z + 1 = 1/ eJ^'^ — 1. In order to compare this scaling 
prediction with the z values obtained directly from the 
finite-size scaling of R, we need to estimate l/eT^*-^. As 
described in Sec. g T'^° is given by T^^ = T/{2ttJ{U")) 
and 1/e = limfc^o l/e(fcjO). However for FTBC we have 



— = lim — -- — r- 

e k^Q (\k, 0) 



1 



> 



e(0,0) 



= 0. 



So for each size L we estimate 1/e by l/e(27r/L, 0). As 
mentioned in the beginning of this section we can also 
include a small correction due to the finite time step 
Ai in the simulations for each size L by replacing T 
by an effective temperature T''^ = (5(0, 0)/(C/"). Fig- 
ure ^ shows Oscaic = z + \ ~ l/eT'^^ — 1 estimated 
in this way as a function of L. When comparing with 
the a values obtained from the finite-size scaling of R, 
we take an average over the relevant L interval. These 
values are shown in Table ||. As is apparent from Ta- 
ble |[ the values of a determined from the size scal- 
ing of R and r agree very well with Oscaio both for the 
RSJ case and the TDGL case. Thus we conclude that 
z = [l/eT'^^) — 2. This conclusion has also been reached 
for the lattice Coulomb gas model with Monte Carlo dy- 
namics.EiJ Furthermore, by invoking the scaling argument 
described above, we infer that the IV exponent should 
be given by a = a^caic = z + 1 = l/eT^^ - iB 

The model given in Ref. || suggests the finite-size scal- 
ing R oc in agreement with our resultsO How- 
ever, according to the reasoning in Ref. H, the scaling ar- 
gument L cx l/j7 leading to the nonlinear IV exponent 
a = ciscaie should break down for small enough currents 
and in this limit one should instead recover a = oahns- 

In the next section we investigate the nonlinear IV 
characteristics more directly by imposing an external cur- 
rent. 



V. NONLINEAR IV CHARACTERISTICS 



In order to obtain the IV characteristics for the 2D XY 
model with RSJ dynamics we use FTBC and Eqs. ( |3^ ) - 
(^). Figure |ll| shows the data obtained from lattice sizes 
L = 4 to 64, where v = V/L is plotted against id = Id/L 
in a log-log plot. As seen from the figure the data are size 
dependent but for L = 64 the data appear to be reason- 
ably size converged except for the smallest currents. The 
data in the figure are for T — 0.80 and the straight line 
is a least-square fit to the L — 64 data in the current in- 
terval 0.03 < id < 0.15 and gives a « 4.7, which is in rea- 
sonable agreement with agcaio = l/eT*-^^ — 1 « 4.5. In the 
following we will investigate the sensitivity of this appar- 
ent agreement to finite size, finite current, and boundary 
conditions. 

One finite current effect is that the exponent a 
refers to a constant Coulomb gas temperature T*-^*^ — 
T /[2-kJ UJJ')]. Since a finite current changes the value 
of (J7"),ll3 fixed temperature (T = const) is not entirely 
equivalent to fixed Coulomb gas temperature (T^"-^ = 
const). In order to convert the data to fixed Coulomb gas 
temperature we have calculated v and T'^^ for T — 0.79 
and 0.80 for fixed external currents, and then by inter- 
polation estimated the voltage value corresponding to 
a fixed T*-^*^. The resulting data for a fixed Coulomb 
gas temperature {T^^ k, 0.17) are shown in the inset 
of Fig. The broken line in the inset is a fit to the 
data and gives a « 4.5. Thus this correction leads to a 
somewhat smaller value of a. 

All previous estimates for the nonlinear IV exponent 
for the RSJ | mO | dq l have been obtained for L = 64 or 
smaller sizes.Bll3l23 The next question we address is how 
much the possible remaining size effects could alter the 
results inferred for L = 64. Figure |l^ shows voltage v 
versus the system size L at the external current id = 0.1 
and T = 0.8 for three different cases. The open squares 
at the top correspond to the usual uniform current in- 
jection method used in Ref. The filled circles corre- 
spond to our FTBC boundary condition and finally the 
open triangles at the bottom correspond to the busbar 
boundary condition used in Ref. It is clear from the 
figure that L = 00 result cannot be estimated by the 
L = 64 for id = 0.1. For smaller id the situation quickly 
gets even worse. Thus this unexpected strong size de- 
pendency clearly makes all earlier results obtained for a 
from IV simulations somewhat questionable .BO'EJ 

As seen in Fig. ^ the uniform current injection ap- 
pears to approach the L — 00 value from above whereas 
the FTBC and the busbar condition appear to approach 
the L = 00 value from below. We have found this to 
be generally true. From this we conclude that L — 256 
is enough to estimate the L = 00 limit for id > 0.1, 
since the data for FTBC and uniform current injection 
are closely the same in this case. The value of a ob- 
tained in this converged current region is about a « 4.1, 
which is somewhat smaller than a « 4.3 obtained from 
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the finite-size scaling of R in the previous section. 

In order to get some further insight, we note that 
the present simulation gives the resistance R = v/id as 
a function of id, as discussed in the previous section, 
for small enough current densities J', XjJ should cor- 
responds to a finite L. Consequently R{c/id), where c is 
a constant, obtained in the present simulations should be 
equivalent to R{L) obtained in the previous section: For 
an appropriate choice of the constant c the data for these 
two simulations should fall on a single curve. Figure |l^ 
illustrates this equivalence, the filled circles are the data 
for R{L) and the open squares are the IV data obtained 
from FTBC with L = 256. The open circles are the aver- 
ages between the L = 256 result for FTBC and uniform 
current injection. When the open circles and squares 
overlap, the L — oo limit has been reached. As seen from 
the figure the two data sets for i? to a good approximation 
fall on a single curve. For large currents R approaches 
the junction resistance r = 1 and for small currents 
R cx {id)""^- The full drawn curve {R = e'-a-'^^Koiu^) 
where Kq is a modified Bessel function) interpolates be- 
tween these two limits [Ko{x) —Inx for small x and 
Ko{0) = 0]. Since the converged IV data are higher up 
on the curve one expects an apparent smaller a than for 
the R{L) data which are lower down on the curve. Our 
conclusion is that the results from the IV simulations 
and the R{L) simulations are consistent with each other 
and with the scaling assumption. 



Scaling collapse 

It is in fact possible to demonstrate the validity of the 
scaling assumption in a more general way: At fixed tem- 
perature R is only a function of L and J^. From the fact 
that R ~ a,t J' = and that the combination J'L 

is dimensionless, one expects that 



R = 



fiJL) 



(42) 



where f{x) is a dimensionless scaling function. The scal- 
ing function f{x) must have the limits /(O) =const since 
R ^ \/L' for J' = 0, and f{x) oc x for large x. The 
latter follows because the L oo limit has to give a 
nonvanishing finite R. This means that the combination 



LR^'" = fiJL) 



(43) 



is only a function of JL. In Fig. |lj we have plotted all 
our simulation data for id < 0.6 as LR^^^ against i^L, 
i.e., the data shown in Fig. ^ together with data for 
L = 128 and 256. Using z as an adjustable parameter, 
we find that all the data collapse onto a single scaling 
curve for z w 3.3. We emphasize that this scaling col- 
lapse involves only one free parameter, z. One also notes 
that the best value for the collapse (obtained by a least- 
square method) is closely the same (z « 3.3 at T = 0.80) 



as was found in the absence of external currents shown 
in Fig. @(a). Furthermore, this zero-id data collapse onto 
a single value for z « 3.3 when plotted as LR^^^ and this 
constant value is given by the broken horizontal line in 
Fig. 14, Thus the data collapse shown in Fig. clearly 



demonstrates that the scaling assumption is valid for all 
the data we have obtained. Since the scaling assumption 
gives a = ascaic = z + 1 — 1/eT^^ — 1, our conclusion 
is that ascaic is indeed the correct IV exponent over a 
broad parameter range. 

The model discussed in Ref. ^ suggests that for small 
enough id the scaling assumption should break down. 
Thus for such small currents the data for large enough 
L should fall above the scaling curve in Fig. |l^. There 
is no sign of any such deviation in our data. However, 
this does not preclude the possibility that such a devia- 
tion could in principle occur for larger sizes and smaller 
currents than we have been able to investigate. 

It is also interesting to note that the scaling function 
f{x) is intimately connected to the finite-size dependence 
of the voltage for FTBC. [See, for example. Fig. |l^ for 
T = 0.8 and id = 0.1 (filled circles).] According to 
Eq. (E^) we have 



f{Ltd) 



Lid 



(44) 



The full drawn curve in Fig. ^ gives u as a function of 
L using Eq. ( ^ ) for id = 0.1 where the scaling function 
f{x) has been obtained by a data smoothing of the data 
in Fig. |l^ The filled circles is a replot of the finite-size 
dependence given as filled circles in Fig. |l2|. As is appar- 
ent from Fig. [ij, the particular shape of the finite-size 
dependence is a direct reflection of the scaling function 

The AHNS predictiono for the nonlinear IV exponent 
differs from the scaling prediction and is instead given by 



1 



OAHNS 



2eTCG 



1. 



The corresponding values are given in Table | and Fig. |T^. 
Our simulations support the scaling prediction. E.g., for 
T — 0.8 and RSJ we find a « 4.3 which is close to the 
scaling prediction Ogcai ~ 4.4 and differs from the AHNS 
prediction cahns ~ 3.7. 



VI. SUMMARY AND COMPARISONS 

The flrst main result of the present investigation is that 
the 2D XY model with RSJ dynamics is well described by 
the MP form for the dynamical response. This appears to 
be generic for 2D vortex fluctuations since the same form 
has-jbeen found for the XY model with TDGL dypam- 
icSjQ the 2D Coulomb gas,|Sdth Langevin dynamical^ as 
well as in experimentsE'HIij However, since the 2D XY 
model with RSJ dynamics is generally accepted as a valid 
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model for a 2D Josephson junction array, the present in- 
vestigation ties the MP form found in the present and 
previous siHiuJations cfoser to the MP form found in ex- 
perimentsE'Ej 

We found the critical exponent z = 2 at the KT 
transition from the finite-size scaling of the resistance R 
using the fluctuating twist boundary condition FTBC, 
both in case of RSJ and TDGL dynamics. Further- 
more, we found the same value of z for RSJ and TDGL 
for all temperatures below the transition using the same 
method. However, we also found that the finite-size scal- 
ing with PBC gave different results. Thus it appears as if 
the finite-size scaling determination of z depends on the 
boundary condition. Our conclusion is that it fails for 
PBC because the characteristic time r is inversely pro- 
portional to the resistance R and for PBC the resistance 
R is identically zero for any finite size. This suggests 
that the proper value of z cannot be determined from 
finite-size scaling with PBC. 

The exponent z determined from the finite-size scal- 
ing with FTBC were found to be the same for RSJ and 
TDGL dynamics and to support the scaling prediction 
in agreement with what was found in 
Ref. |l^ for the 2D lattice Coulomb gas with Monte Carlo 
dynamics. We also explicitly showed that the exponent 
z determined from the finite-size scaling of R is related 
directly to a diverging relaxation time. Thus our con- 
clusion is that z is larger than 2 below the KT transi- 
tion. This_cpsult is in agreement with thjajnodel discussed 
in Ref. Using a scaling argument, El we related the 
finite-size scaling of R to the nonlinear IV characteris- 
tics by noting that the current density plays the role 
of 1/L leading to V (x 1°" with a — z + 1. Consequently, 
provided the scaling argument is valid,|-our simulations 
support the prediction a — — 1.E3 

We also calculated the IV exponent a directly from 
the voltage ^ as a function of current /. Here we found 
that the results were strongly size dependent. This large 
size dependence we found for standard current injection 
boundary, FTBC, and the "busbar" boundary condition 
introduced in Ref. ||. For our largest lattice sizes 256 x 256 
a size-converged result could only be estimated for cur- 
rents which seemed to be outside the true scaling regime 
V (X 1°". However, by using the relation L (x \/J valid 
for small enough J we showed that the data for the re- 
sistance simulation R{L) and the IV simulations R{c/ J) 
can be made to fall on a single curve for an appropriate 



choice of the constant c. This agreement suggests that 
our IV simulations and our R{L) simulations are con- 
sistent with each other and with the scaling assumption. 
We concluded that it is difficult to obtain the nonlinear 
IV exponent a directly from the V{I) data in case of the 
2D XY model with RSJ dynamics. This is because resis- 
tance ratios R{I)/r < 10^^ (r is the junction resistance) 
seem to be needed. This in turn implies such small cur- 
rents that lattice sizes considerably larger than 256 x 256 
are required to avoid the finite-size effects. However, m 
case of the 2D Coulomb gas with Langevin dynamicst2l 
it has been possible to converge the simulations closer to 
where the true scaling V oc I'^ appears to be valid and 
in these cases the scaling exponent a — — 1 was 

deduced from the V{I) data. 

Finally, we showed that all our IV data and our R{L) 
data for a fixed temperature collapse onto a single scal- 
ing curve f(x = Lid). This data collapse demonstrates 
that the scaling argument is indeed valid over a broad 
parameter range and thus confirms that the nonlinear 
IV exponent is given by flgcaic = l/eT^*^ — 1 over the 
parameter range covered by our data. This does not pre- 
clude the possibility that, for smaller currents and larger 
sizes than we have been able to converge, there might 
be a deviation from the scaling curve given in Fig. |l^ as 
suggested by the model in Ref. ^. However, there is no 
sign of any deviation from the scaling curve in our data 
for the RSJ model. 

In short, the present simulations of the 2D XY model 
with RSJ dynamics confirm the picture that 2D vortex 
fiuctuations has an anomalous kind of dynamics. The 
characteristic features of this dynamics are presumably 
linked to the logarithmic vortex interaction. However, 
a firmer theoretical understanding of the characteristic 
features, which have been encountered in numerous sim- 
ulations, as well as in experiments, is still lacking and is 
a challenge for future research. 
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APPENDIX: LINEAR RESPONSE 



A total current ix{Y,t) which varies slowly in time compared to the thermal fluctuations gives rise to an average 
nonvanishing phase difference q{r,t) — {S/ xO{r,t)) . Thus Eqs. and (|lj) together with the chain rule gives 



P(r 



,t~t') 



-J 



cPr"dt 



„ a(C/'[V,0(r,t)]) 



ag(r",t") 



9(V,0(r",t")) 



+ 6{r-r')S{t-t'), 



(Al) 
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where |o denote that the resulting averages should be the equilibrium ones. Let us introduce the notation 

d{U'[VJ{v,t)]) 



Q{v-r",t-t") = J 



dq{v",t") 



then the Fourier transform of Eq. ( Al) is just 

ia;P(k,u;) = -Q(k, w)P(k, w) + 1 

so that 

1 



(A2) 



P(k,w) 



iuj + (5(k, w) 



We note that 



Q(r-r',t-i') = J 



d{u'[vj{v,t)]) 



dq(r\ t') 



J{U"[V0{r, t)])S{r - r')Sit ~t') + J 



d{U'[V,e{r,t)]) 



dq(r',t') 



Here the last term is for t ^ t' and r 7^ r' so that the disturbance q{r',t') = (Va;6'(r', i')) couples linearly to 
JU'[Vxd{r' ,t')] in the XY Hamiltonian. Consequently, the corresponding correlation function is 

-j2([/'[v,0(r,t)][/'[V,0(r',t')]) 
and by the fluctuation-dissipation theorem we have 

Q{r,t) = ^^({/'[V,0(r,i)][/'[V,0(O,O)]) + J(C/"[V,0(O,O)])<5(r)^(i) 

for t > and otherwise. Next we note that a space Fourier transform of the correlation function U'[\/xd{r,t)] 
U'[\7 x^iic' ,t')]) gives the correlation function G(k, t) defined in connection with Eq. (||) so that 



TJo 



3 - 

dte-*-*-G(k,i) 



T 



T 



e(k,w)' 



(A3) 



where po = J{U") and the result is obtained by partial integration and comparison with Eqs. (^-(|6|). 
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TABLE I. Comparison between the exponent a = z + 1 
obtained from the R{L) simulations and the predicted values 
tSscaic and oahns for RSJ and TDGL dynamics. The values of 
Sscaic and flAHNS are obtained from the averages over L = 10, 
12, and 16 (see Fig. |l^ for RSJ case). The exponent a = z + 1 
is obtained from R{L) ~ in Fig. |^ and is found to be 
consistent with 2: in r ~ in Fig. ^. The numbers in paren- 
theses represent the statistical errors of the last digits. It is 
clearly shown that the exponent a measured by direct calcu- 
lation of resistance from Eq. (^l|) is much closer to Oscaic than 
to flAHNS for both RSJ and TDGL dynamics. 
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FIG. 1. The dynamical response function l/e{0,uj) of the 
2D XY model with RSJ dynamics at T = 0.85 for a 64 x 64 
lattice with periodic boundary conditions. [The frequency u 
is in units of 2eric/h (see text).] The filled squares and cir- 
cles correspond to the real part and the absolute value of the 
imaginary part of the dynamical response function, respec- 
tively. The full curve is obtained by fitting to the real part 
of the MP form response function in Eq. (W and the broken 
curve is the imaginary part Eq. using the same values 
of the fitting parameters as for the full curve. The vertical 
broken line corresponds to the u for which the peak ratio 
|Im[l/e(0,cj)]|/Re[l/e(0,t^) - 1/^(0,0)] is 2/tt. At this lu the 
absolute value of the imaginary part should, accordingly to 
the MP form, have a maximum. 




f 

FIG. 2. The time correlation function ln[fG(0, t)] versus 
time t at T = 0.85 for various system sizes [L = 6, 8, 10, 12, 16, 
and 64 from bottom to top] in case of (a) RSJ and (b) TDGL 
dynamics. The full curves have been obtained by fitting to 
the modified-MP (MMP) form Eq. (|3|). The figure shows 
that that the relaxation time tq in the MMP form diverges 
as the system size is increased and that (5(0, t) on 1/t for large 
t in the thermodynamic limit. 
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RSJ 
TDGL 
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FIG. 3. The relaxation time tq for the time correlation 
function 6(0, t) a.t T = 0.85 for RSJ (squares) and TDGL 
(circles) dynamics. The data points have been obtained from 
least-square fits to the MMP form G^^p given by Eq. (^ as 
shown in Fig. |^. As the system size L is increased tq diverges. 
However, the exponent z defined by tg ~ appears to have 
different values for the two types of dynamics. The lines are 
obtained from least-square fits using data points for L = 8, 
10, 12, and 16 in the RSJ case and L = 6, 8, 10, and 12 in the 
TDGL case, giving z fa 1.6 and z ^ 2.0 for RSJ and TDGL, 
respectively. 
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FIG. 4. The dynamical response function l/e(k, oj) at 
k = k,„in = (0,27r/L) for a L = 64 array at T = 0.85 for 
RSJ dynamics with FTBC (cu is in units of 2eric/h). The 
filled squares and circles correspond to the real part and 
(the absolute value of) the imaginary part of l/e(k,Lj). The 
full curve is obtained from a least-square fit to Eq. (H) with 
two free parameters ljq and e and the broken curve is ob- 
tained from Eq. (^ using these parameter values. The ver- 
tical broken line is given by the condition that the peak ra- 
tio \Ira[l/e{k,Lj)]\/Re[l/e{k,uj) - l/e(A:,0)] = 2/7r and at this 
value of Lu, the absolute value of the imaginary part should, 
according to the MP form, have a maximum. 
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FIG. 5. The dynamical response function l/e:(k, w) at fi- 
nite k for a L = 64 array at T = 0.85 for RSJ dynamics with 
FTBC {(jO is in units of 2eric/Ti). The real and imaginary parts 
are obtained using the wave vectors k = (0, fcy = 2-Kny/L) 
with Uy = 1, 2, A, 6, 8, and 10 (from top to bottom). The 
imaginary part depends very little on the value of k in the 
frequency interval around the maximum, in contrast to the 
real part which increases with decreasing k. 




FIG. 6. (a) The time correlation function ln[iG(k, t)] ver- 
sus time t eXT = 0.85 for RSJ dynamics with FTBC. The 
wave vectors are k = (0, ky = 2'Kny/L) with ny = 1, 2, 4, 6, 
8, and 10 (from top to bottom) and the array size is L = 64. 
As k ^ 0, G(k,t) approaches G(k,t) 1/t for large values 
oft. At nonzero value of k, G{k,t) exhibits exponential decay 
G{k,t) ~ exp[—t/TG{k)]/t in the long-time limit. The broken 
lines are plotted with the Taik) values corresponding to the 
straight line in (b), where we show rcik) versus k^ . In (b), 
the squares have been obtained from the least-square fit of 
G{k, i) to the exponential decay form, and the full straight 
line is the result of the least-square fit to tg {k) oc k^ . 



17 



10" 



R 



7=0.80 □ 
7=0.85 o 
7=0.90 




LA{t) 




20000 



40000 



60000 80000 100000 



10-^ 


(b) 

A 


7=0.80 □ 
7=0.85 
7=0.90 


10-2 

R 






10-3 








6 8 


10 12 16 



FIG. 8. Time evolution of the variable LA{t) at T = 0.85 
as a function of time t for RSJ dynamics and system sizes 
L — 6, 8, 10, 12, and 16 (from top to bottom). The curves 
are shifted in the vertical direction. As seen LA{t) sometimes 
makes discrete jumps of size 2-k (the unit of the vertical axis 
is 27r). The characteristic time r in Fig. ^ is related to the 
average time between the 2-k jumps. 



FIG. 7. Resistance 7? versus system size L for (a) RSJ and 
(b) TDGL dynamics obtained from Eq. (^. The full lines 
are obtained by fitting to the scaling form i? ~ and from 
these fits the values of z are determined to be z = 3.3(1), 
2.7(1), and 2.0(1) at T = 0.80, 0.85, and 0.90 for the RSJ 
case, and z = 3.3(1), 2.8(1), and 2.1(1) at T = 0.80, 0.85, and 
0.90 for TDGL. 
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FIG. 9. The relaxation time r obtained directly from the 
time scale of the 27r jumps of LA{t). The obtained values of 
T are plotted against the system size L for (a) RSJ and (b) 
TDGL dynamics (see Fig. fel). The full lines represent t ^ 
with the z values taken from Fig. |^. The figure illustrates that 
z determined from the scaling of the resistance R is indeed 
associated with a divergent characteristic time. 
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FIG. 10. Predictions of the IV exponent for the RSJ model 
at T = 0.8 as a function of the system size L. The open 
squares are obtained from Oacaie = — 1 for FTBC 

whereas the open circles represent oahns = l/2?riif + 1 for 
FTBC. 
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FIG. 11. The current-voltage (JV^ characteristics at 
T = 0.8 for the fluctuating twist boundary condition at var- 
ious system sizes. The full straight line is obtained from the 
least-square fit in the interval 0.03 < id < 0.15 for L = 64 
which gives a « 4.7. The linear region with IV exponent 
a = \, seen for the smaller sizes and small currents (the dot- 
ted straight line has the slope a = 1), disappears as the system 
size is increased. Inset: IV curve for L = 64 at fixed Coulomb 
gas temperature t'^'^ « 0.17, corresponding to T = 0.80 with 
no external currents. The broken line is obtained from the 
least-square fit in the interval 0.07 < id < 0.15, giving a ~ 4.5. 
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FIG. 12. Voltage v versus system size L at the current 
id = 0.1 for T = 0.8. The empty squares are for the uni- 
form current injection with periodic boundary conditions in 
the direction perpendicular to the current. The empty trian- 
gles are obtained with the critical current ic = 10 for vertical 
junctions on the boundaries, which is very similar to the bus- 
bar boundary. The filled circles are for FTBC introduced in 
Sec. III. As the system size is increased, the voltages for all 
three methods are shown to converge towards the same value 
in the L — oo limit. However, the uniform current injection 
approaches the L = oo limit from above whereas the FTBC 
and busbar condition approach from below. The lines are 
guides to the eye. 



FIG. 13. The resistance R{c/id) = v/id at T = 0.80 (open 
squares correspond to the L = 256 data for FTBC and open 
circles to the average between the L = 256 data for FTBC 
and the uniform current injection) is compared to the resis- 
tances R{L) at T — 0.80 [filled circles, the same data as in 
Fig. ^(a)]. Choosing the constant c ~ 0.70 makes the two 
data sets collapse onto a single curve. The full drawn curve 
interpolates between the limits R = 1 for large currents and 
R cc {id)'^~^ for small currents (the explicit form of the in- 
terpolation curve is ef^-^)^'"''/^) with a = 4.3 and b = 1.42). 
The figure suggests that the two ways of calculating R are 
consistent and that the R{c/id) data are not quite in the 
asymptotic small-current regime. 
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FIG. 14. Demonstration of the validity of the scaling as- 
sumption. The IV data for the RSJ model with T = 0.8 and 
id < 0.6 are plotted as LR^^^ against Lid- For z « 3.3 all the 
data for the various L and id collapse onto a single scaling 
function f{x = Lid)- The horizontal broken line corresponds 
to the constant value for LR{LY^^ obtained for id = for the 
same value of z [see Fig. ^(a)]. The straight line corresponds 
to the linear behavior f{x) ~ x for large x. 
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FIG. 15. The relation between the finite-size dependence 
of the voltage v and the scaling function f(x — Lid)- The 
full drawn curve is the function v = id'^^[f{x)/xY where f(x) 
has been obtained by a data smoothing of the data in Fig. [14. 
The filled circles are the finite-size data for w at T = 0.8, the 
same data as the filled circles in Fig. [13. 
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